城市 大 气 环 境 的 高 保 真 度数 值 模拟 
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随 着 城市 人 口 和 工业 活动 的 增加 ， 我 国 多 数 一 
二 线 城市 在 建 及 规划 建筑 呈 高 层 化 、 密 集 化 趋势 。 
其 对 大 气 环境 最 直接 的 影响 是 大 幅 降低 城市 内 空气 
流通 速度 。 空 气 停滞 形成 的 空气 穹 隆 促使 近 地 逆 温 
层 的 形成 ， 将 各 种 污染 源 排放 的 污染 物 笼 置 ， 加 剧 
城市 的 污染 ， 形 成 高 浓度 的 大 气 污染 [1]。 此 外 ， 城 
市 大 气 环境 特 性 的 变化 ， 如 因 建 筑 物 阻挡 产生 的 涡 


2) ”为 城市 布局 及 规划 提供 依据 

为 应 对 各 大 城市 日 趋 严峻 的 大 气 污染 状况 ， 除 
了 从 源头 上 减轻 污染 源 的 措施 外 ， 城 市 风 道 作为 治 
理 大 气 污染 的 手段 之 一 成 为 研究 的 热点 领域 。 大 气 
污染 物 在 静 稳 无 风 、 大 气 扩散 条 件 差 的 不 利 气象 条 
件 下 容易 堆积 ， 因 此 把 郊外 的 风 引 进 主 城区 ， 将 性 
等 污染 物 吹 走 成 为 备 选 除 考 方 式 之 一 。 城 市 风 道 即 
城市 通风 廊 道 ， 通 过 以 市 区 中 的 绿地 、 水 面 及 道路 
等 形式 减少 城市 对 风 的 阻挡 ， 达 到 促进 城区 内 外 热 
交换 和 污染 物 扩 散 的 作用 ， 从 而 缓解 热岛 效应 和 空 


旋 大 大 增多 、 气 流 的 素 流 性 增强 等 ， 会 在 城市 街道 
中 产生 “城市 急流 ”或 产生 气流 死 区 ， 对 市 区 街道 、 
小 区 等 小 面积 区 域 的 空气 污染 物 扩散 带 来 不 利 影响 
甚至 造成 严重 污染 问题 [2]， 同 时 也 会 使 城市 大 气 环 
境 数值 计算 变 得 困难 。 


2 需求 


在 污染 预报 等 环境 监测 领域 ， 基 于 背景 流 场 的 


气 污染 。 北 京 市 政府 从 2014 年 开始 分 析 、 论 证 与 规 
划 通 风 廊 道 ， 以 增强 通风 潜力 、 绥 解 热 岛 效 应 。 据 
媒体 报道 ， 北 京 将 形成 5 条 宽度 500 米 以 上 的 一 级 通 
风 底 道 ， 多 条 宽度 80 米 以 上 的 二 级 通风 底 道 ， 未 来 
形成 通风 亡 道 网 络 系统 。 相 关 研 究 建议 ， 对 主 通风 
亡 道 区 域 严格 规划 控制 ， 包 括 控制 建设 高 度 和 密度 
等 ， 同 时 打通 障碍 点 。 除 北京 外 ， 上 海 、 杭 州 、 武 
汉 等 多 个 城市 也 开始 着 手 规划 城市 风 道 ， 以 应 对 越 
来 越 严重 的 大 气 污染 问题 。 在 城市 规划 过 程 中 ， 模 


模拟 与 相关 部 门 的 气象 数据 ， 求 解 大 气 污 染 物 输 运 
及 扩散 问题 获取 的 数值 结果 可 作为 空气 污染 预报 的 
重要 参考 资料 ， 可 用 以 评估 大 气 排放 的 污染 物 对 城 
市 区 域 环 境 造成 的 影响 。 近 年 来 ， 在 城市 突 发 性 污 
染 及 危险 事件 分 析 及 预警 、 城 市 规划 等 方面 ， 大 气 
环境 的 高 保 真 度数 值 模拟 应 用 越 来 越 广 。 

1) 城区 突 发 污染 事件 预报 与 预警 

快速 有 效 地 对 城市 大 气 环境 进行 模拟 ， 并 以 此 
安排 预报 预警 措施 具有 重大 的 社会 和 经 济 意义 。 以 
深圳 为 例 ， 该 市 及 周边 地 区 突 发 污染 事件 如 油库 、 
化 学 品 仓库 等 高 危 品 爆炸 ， 工 业 厂 区 发 生火 灾 等 时 
有 发 生 。 如 与 深圳 市 中 心 城区 域 接壤 的 香港 打鼓 岭 
垃圾 回收 场 分 别 与 2013 年 11 月 9 日 、2015 年 9 月 11 日 、 
2016 年 3 月 2 日 三 次 发 生火 警 ， 大 量 夹杂 有 毒气 体 的 
黑 烟 球 向 深圳 方向 , 对 城区 的 空气 造成 严重 的 污染 ， 
严重 威胁 市 民 的 身体 健康 。 基 于 高 保 真 度 的 大 气 背 
景 流 场 进行 污染 物 实时 扩散 及 输 运 的 数值 模拟 ， 可 
及 时 对 污染 物 在 市 区 各 区 域 的 污染 浓度 做 出 计算 ， 
获取 ， 为 进一步 的 评估 、 预 警 及 其 他 调控 治理 措施 
提供 依据 。 


拟 大 气 环境 以 及 污染 扩散 状况 可 以 为 政府 部 门 提供 
各 种 有 价值 的 参考 信息 ， 合 理 规 划 建 筑 物 的 密度 和 
高 度 以 利于 空气 污染 物 、 城 市 内 部 热量 和 废弃 物 朴 
散 ， 在 为 城区 保留 安全 绿色 空间 和 保持 生态 格局 前 
提 下 优化 城市 布局 。 

城区 空气 质量 与 人 民生 活 质量 恩 息 相关 ， 市 中 
心 区 人 口 分 布 又 非常 密集 ， 因 此 及 时 、 准 确 地 对 城 
区 大 气 环境 及 基于 其 基础 上 的 污染 预报 与 预警 ， 以 
及 以 大 气 污染 消除 与 控制 为 目的 的 合理 的 城市 规划 
对 于 环境 保护 及 市 民 健康 意义 重大 ， 世 界 各 国 也 都 
在 积极 地 开展 城市 大 气 环境 的 研究 。 


3. 研 究 现状 与 挑战 


D ”研究 现状 与 趋势 

运用 数值 模拟 技术 对 城市 的 大 气 环 境 进 行 模 
拟 已 经 成 为 相关 领域 的 研究 热点 .在 模拟 模式 方面 ， 
近 几 十 年 来 通过 国内 外 相关 学 者 的 大 量 研究 ， 以 梯 
度 输 送 、 满 流 统计 和 满 流 相似 扩散 三 种 理论 为 基础 ， 
研究 者 已 经 建立 和 发 展 了 多 种 大 气 环境 模拟 模式 ， 


主要 有 高 斯 模式 及 其 变形 模式 、 统 计 模 式 、 大 气 扩 
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呈 几 何 多 尺度 , 离散 后 非 线性 系统 病态 且 规 模 巨 大 。 


散 相 似 模式 等 。 相 应 的 大 气 环境 模型 软件 也 已 经 开 
发 出 数 百 种 ， 如 美国 EPA 的 Model-3 与 
ISC-AERMOD [3]、 英 国 CERC 的 ADMS [4] 等 。 不 
同 的 扩散 模式 及 与 其 相应 的 模型 软件 都 有 各 自 不 同 
的 地 貌 、 气 候 、 污 染 源 等 考虑 因素 ， 因 而 也 就 有 各 
自 不 同 的 适用 范围 。 城 市 大 气 环境 模拟 从 模拟 尺度 
方面 可 分 为 如 下 几 类 [5]: 第 一 类 是 城市 冠 层 模式 模 
拟 ， 分 辨 率 数 百 至 数 千 米 ， 预 报时 间 为 数 小 时 。 第 
二 类 是 单个 到 多 个 建筑 物 尺度 的 模拟 ， 典 型 网 格 尺 
度 为 数 米 至 数 百 米 ， 采 用 CFD 模 式 ， 求 解 时 间 约 为 
数 分 钟 至 数 小 时 。 第 三 类 为 小 于 建筑 物 尺 度 的 模拟 ， 
典型 网 格 尺度 小 于 一 米 ， 同 样 采用 CFD 模 式 ， 求 解 
时 间 为 数秒 至 数 分 钟 。 

由 于 计算 机 硬件 条 件 限制 , 过 去 大 部 分 的 研究 
往 需 要 首先 确定 数值 模拟 的 尺度 问题 ， 根 据 问题 规 
模 选 择 不 同 的 模拟 方法 ， 如 中 尺度 模拟 或 小 尺度 模 
拟 。 其 中 ， 中 尺度 问题 虽然 可 极 大 减 小 计算 量 ， 牺 
牲 仿真 的 精度 来 换取 计算 的 可 实现 性 ， 但 还 是 简化 
城市 下 垫 面 建筑 模型 ， 结 果 相对 于 特定 建筑 单元 是 
定性 而 非 定量 结果 ， 因 此 对 大 气 环 境 模 拟 的 精度 较 
差 。 采 用 CFD 方 法 的 小 尺度 污染 数值 模拟 同样 受制 
于 计算 机 硬件 的 水 平 ， 只 能 获取 某 栋 或 数 栋 建筑 物 
周边 大 气 环境 情况 ,无 法 与 该 区 域 整体 环境 相 结合 ， 
因此 往往 准确 度 不 高 ， 仅 在 室内 等 较 封 闭 空间 内 可 
以 获取 可 置信 的 数值 结果 。 近 年 来 ， 城 市 大 气 环境 
中 小 尺度 模拟 发 展 趋势 主要 概括 为 以 下 三 点 : 首先 
是 模拟 方案 分 辨 率 更 高 ， 能 更 精确 表达 城市 地 面 的 
复杂 性 和 相互 关系 ， 中 尺度 与 小 尺度 问题 紧密 结 
合 ， 模 拟 更 细致 、 更 新 颖 ;第 二 是 可 在 更 大 、 更 复 
杂 区 域 实施 模拟 ; 第 三 结果 评价 更 严格 、 更 细致 ， 
可 信和 度 提 高 [6]。 推 动 城市 大 气 环境 中 小 尺度 模拟 克 
合 相 关 研 究 发 展 的 两 个 重要 因素 是 城市 人 口 持续 增 
长 造成 的 环境 问题 愈加 严重 以 及 计算 机 硬件 技术 的 
发 展 使 数值 计算 能 力 快速 增长 [7]。 

2) ”挑战 性 问题 

城市 的 典型 尺度 一 般 认 为 约 20~40km， 城 市 下 
热 面 建筑 群 即 边界 层 按照 功能 与 用 途 可 不 均匀 地 划 
分 成 许多 区 域 ， 包 括 市 中 心 区 、 居 民 区 、 公 共 活 动 
区 、 新 建 区 、 工 业 区 、 不 同 街 道 交 通 区 等 。 一 方面 ， 


高 分 辨 率 的 数值 结果 往往 需要 高 精度 的 流动 模型 以 
及 高 分 辨 率 的 时 空 网 格 ， 这 些 因素 均 带 来 了 庞大 计 
算 量 ， 求 解 时 需要 超 强 的 计算 机 硬件 与 与 之 匹配 的 
高 性 能 数值 算法 。 

数值 模拟 技术 的 发 展 始终 与 计算 机 硬件 发 展 水 
平 紧密 相关 。 对 城市 大 气 环境 进行 模拟 ， 实 质 上 是 
以 气象 预报 和 监测 信息 为 初 边 值 条 件 ， 对 空气 运动 
的 数学 方程 一般 是 三 维 非 线性 偏 微 分 方程 组 ) 进 行 
计算 机 求解 。 虽 然 此 类 数值 方法 已 经 应 用 多 年 ， 并 
已 经 是 高 性 能 计算 最 重要 的 应 用 领域 之 一 ， 但 由 于 
传统 的 大 气 环 境 模拟 并 行 算法 及 其 程序 优化 技术 不 
能 适应 目前 集群 处 理 器 数目 高 速 增长 的 需要 ， 以 至 
该 领域 几乎 所 有 程序 无 法 扩展 到 数 干 核 ， 韦 论 千 万 
亿 次 的 数 万 核 乃 至 亿 亿 次 机 的 数 十 万 、 数 百 万 核 。 
上 述 原因 致使 目前 城市 大 气 环境 数值 模拟 相关 应 用 
缺乏 高 扩展 并 行 算 法 的 支撑 ， 直 接 导 致 数值 模拟 求 
解 分 辩 率 大 于 百 米 且 在 万 亿 次 泽 点 运算 量 级 的 计算 
平台 上 展开 ， 严 重 影响 数值 模拟 的 准确 率 。 近 年 来 
天 津 、 深 圳 、 广 州 等 地 相继 建立 了 国家 超 算 中 心 ， 
所 采用 的 超级 计算 集群 普遍 采用 新 型 架构 体系 ， 这 
些 都 对 当前 的 通用 算法 技术 提出 了 严峻 的 挑战 。 而 
研究 适用 于 求解 大 规模 非 线性 稀 玻 方程 组 的 高 可 扩 
展 并 行 算法 与 计算 机 硬件 资源 匹配 已 经 成 为 工程 计 
算 领 域 非常 重要 的 研究 方向 。 因 此 ， 若 获取 及 时 、 
准确 的 污染 物 输 运 扩散 数值 结果 ， 必 须发 展 一 套 能 
够 充分 发 挥 当前 计算 资源 ， 尤 其 是 大 规模 并 行 计 算 
平台 的 高 效能 计算 优化 理论 和 技术 ， 使 得 预报 程序 
具有 高 可 扩展 性 ， 从 而 确保 计算 的 高 速 、 高 保 真 度 。 


4 ”关键 技术 


4. 1 流体 控制 方程 及 其 有 限 元 离间 

为 获取 高 保 真 度 的 城市 大 气流 场 信息 ， 首 先 
对 大 气流 动 以 及 所 研究 的 市 区 建筑 群 进行 建 模 。 
次 , 由 于 流体 的 连续 介质 特性 以 及 边界 条 件 的 影响 ， 
往往 需要 将 建筑 群 置 入 某 个 封闭 的 计算 区 域名 ， 然 
后 针对 整个 2 进行 分 析 求 解 。@ 的 选取 与 计算 问题 
有 关 ， 电 必须 足够 大 以 尽量 避免 区 域 边界 对 建筑 群 
附近 流 场 的 影响 , 同时 还 要 兼顾 整体 计算 量 的 大 小 。 


并 强 


整个 城市 的 大 尺度 性 与 其 中 各 城市 单元 的 几何 多 尺 


最 后 ， 市 区 大 气流 动 可 视 为 不 压缩 Newton 流 ， 可 


度 性 给 高 保 真 度 模 拟 带 来 巨大 的 计算 量 ; 男 一 方面 ， 
边界 层 结构 与 特性 的 非 均匀 性 带 来 了 气象 及 流动 参 
数 难 以 确定 、 动 量 、 热 量 输 运 过 程 复杂 多 变 等 难题 
[8]。 从 计算 角度 考虑 ， 城 市 大 气 环境 模拟 是 一 个 典 
型 的 高 雷诺 数 时 间 相 关 三 维 问题 ， 求 解 区 域 复杂 且 


采用 三 维 非 定常 不 可 压 Navier-Stokes 方程 描述 : 


Ry 
> x (|) 


其 中 ?、 了 及 《分别 为 空气 密度 、 压 力 和 动力 粘度 
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系数 ， CCLVTVB 为 速度 向 量 ， 字 7432 为 柯 西 


压力 张 量 ， 其 中 ，I 为 特征 张 量 ，< 人 二 see 。 
该 连续 方程 及 动量 守恒 方程 由 初始 条 件 : 


本 = ECAP 
以 及 边界 条 件 : 
mm 一 pi Lz - 
恒生 一 和 E》 < 二、 总 
RE < 四 [ER 


约束 ， 其 中 Tm 为 速度 进口 边界 ，Tw 为 城市 建筑 群 
落 及 计算 区 域 边界 ， 设 为 无 滑 移 壁 面 边 界 ，Tow 为 
采用 远 场 边 界 条 件 即 自由 出 流 边 界 ， 为 进口 边界 
条 件 给 定 的 风速 。 在 计算 实际 问题 时 ， 边 界 条 件 由 
气象 及 环保 部 门 提供 的 实时 监测 数据 选 定 。 城 市 建 
筑 群 几何 模型 可 根据 GIS 技术 、 点 云 建 模 技术 以 及 
其 他 建 模 方 法 获取 。 如 图 一 所 示 ， 实 际 模拟 中 ， 处 
理 几何 模型 相当 繁琐 ， 为 获取 高 保 真 度数 值 结果 ， 
往往 需要 将 地 形 地 势 考虑 在 内 ， 其 工作 量 可 占 整 体 


模拟 时 间 的 60% 以 上 。 


妇 一 深圳 会 展 中 心 建筑 群 分 布 及 建 模 示意 图 

为 了 获取 流体 流动 数值 解 ， 数 值 模拟 技术 是 通 
过 将 控制 偏 微分 方程 按照 茶 种 离散 格式 离散 ， 然 后 
转化 为 各 个 计算 节点 上 的 代数 方程 组 ， 最 后 求解 得 
到 各 节点 上 的 值 来 获取 整体 流 场 的 近似 解 ， 因 此 离 
散 化 与 代数 化 是 数值 计算 的 前 提 。 对 于 本 文 的 瞬 态 
问题 , 偏 微 分 方程 的 离散 在 空间 和 时 间 上 都 有 体现 ， 
选取 合适 的 时 空 离散 格式 对 后 续 求 解 算法 有 重要 影 
响 。 


空间 离散 格式 方面 ， 考 虑 到 高 阶 的 有 限 元 在 离 
散 化 后 所 形成 的 矩阵 非 零 元 素 过 多 从 而 影响 算法 整 
体 的 并 行 可 扩展 性 ， 结 合并 行 计算 的 特点 ， 我 们 采 
用 低 阶 的 #-# 有 限 元 格式 对 空间 进行 离散 。 

通过 对 计算 区 域 Q 进行 剖 分 得 到 非 结 构 四 面 
体 网 格 下 ={K}， 定 义 测试 函数 及 权 函 数 空间 : 
U={u(,t) |u(,t) EeE [HION, u(t)=g on Dict 
w= {u(t) | u(t) EHOF, u(t)=0 on TineU Twan} 
P={p(,t) |7(,t) €E LO)}, 

Navier-Stokes 方程 (1) 的 弱 形 式 可 描述 为 寻求 
ueU 及 pe P 使 得 所 有 $eWU，w e P 满 足 ， 


ph Bd +nufo vu: ved +p holu: Vu: Badan 
— Jopv :Bd + ov aedo = hf Gdn 


有 限 维 试 函数 及 全 函数 空间 可 以 表示 为 : 

Ur = {u(t) | t) = Di Ba,t), w(t) =g on Tour 

Ut = {u(t) | ut(,t) = De Bh t), wt)=0 on TiwaUTwan} 
Pr = {p(t) p(t) = Di php t)}, 

其 中 Nu 忆 为 速度 与 压力 的 总 节点 数 , WE 信 ， PY ER 

为 节点 处 速度 与 压力 值 。 由 于 ?43 元 不 满足 

Ladyzenskaja-Babuska-Brezzi (LBB) 条 件 , 因此 ， 

通过 添加 适当 的 稳定 项 “， 空 间 离散 有 限 元 系统 可 

表示 为 ; 

ph .Bd tp vu :ved + pf Vu Gdn 

一 Jo PV BidN 十 人 JCY .Un)eondo 十 2 Ke (V .th TV . B") 

号 ( 竖 + (Ww V+ vp, Tian VB Vp)) 

= hf BrdQ + Der (Ff, mw VE + Vp!)) 


时 间 项 方面 ， 本 文采 用 了 全 隐 时 间 格 式 离散 进 
行 离散 。 全 隐 格 式 由 于 减轻 甚至 消除 了 CFL 条 件 对 
时 间 步 长 的 限制 ， 在 进行 大 规模 数值 计算 时 具有 一 
定 优势 。 有 具体 地 ， 针 对 上 述 不 可 压 Navier-Stokes 
方程 ， 选 取 隐 式 向 后 Euler 有 限 差 分 格式 对 时 间 项 
离散 ， 即 将 控制 方程 (1) 离散 为 如 下 非 线 性 代数 方 
程 : 


rn—rn-l1 ,nN 
| (2) 


其 中 为 时 间 步 长 ，S(tz 为 在 第 "时 间 步 经 空间 
离散 后 形成 的 半 离 散 系 统 。 至 此 ， 问 题 化 解 为 在 每 
个 时 间 步 求解 非 线性 方程 组 : 
F"(2")=0 (3) 

4. 2 并 行 Newton-Krylov-Schwarz 算法 

包括 城市 大 气 环境 流 场 计算 等 问题 在 内 ， 计 算 
流体 力学 技术 的 核心 是 求解 流体 控制 方程 经 有 限 元 
等 方法 离散 后 所 形成 的 稀 玻 方程 组 。 计 算 精 度 要 求 
提高 、 城 区 下 执 面 建筑 建 模 精细 化 、 多 尺度 以 及 大 
雷诺 数 对 流 占 优等 特点 ， 使 得 该 稀疏 方程 组 的 规模 
十 分 巨大 且 伴 随 着 强烈 的 非 线 性 ， 求 解 异 常 困 难 。 
具体 到 城市 大 气 环 境 模拟 问题 ， 在 每 个 时 间 步 内 ， 
经 有 限 元 离散 后 形成 的 系统 (3) 规模 十 分 巨大 且 由 
于 对 流 占 优 而 具有 很 强 的 非 线 性 。 求 解 此 类 大 规模 
的 病态 问题 时 ， 一 般 均 需 借 助 大 型 并 行 计算 机 及 相 
应 的 可 扩展 性 算法 进行 求解 。 

目前 非 线性 方程 组 线性 化 方法 中 最 常用 且 有 效 
的 是 Newton 型 方法 。 其 中 非 精 确 Newton 方法 虽然 
需要 求解 原始 的 Jacobian 矩阵 , 但 是 却 不 需要 精确 
求解 ， 可 通过 求解 精度 来 达到 计算 量 与 收敛 速度 的 
平衡 。 在 线性 求解 器 方面 ， 由 于 直接 法 求解 大 规模 
线性 方程 组 对 计算 机 内 存 及 速度 要 求 过 高 而 不 现 
实 ， 因 此 目前 使 用 友 代 法 求解 大 规模 非 线性 系统 几 
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乎 成 为 唯一 选择 .本 文 基于 Krylov 子 空间 的 迭代 方 
法 构造 线性 求解 器 。 上 述 两 种 方法 结合 头 
Newton-Krylov 迭代 法 ， 该 方法 的 核心 是 在 每 个 非 
线性 迭代 步 通 过 采用 Krylov 子 空间 迭代 法 求解 一 
个 雅 可 比 线性 方程 组 从 而 获得 Newton 增 量 。 
无 论 在 线性 、 非 线性 迭代 法 中 ， 预 条 件 子 都 起 
到 加 速 收敛 的 作用 ， 在 数值 算法 构造 中 占据 重要 地 
位 。 其 作用 是 以 较 低 的 额外 代价 换取 迭代 次 数 的 大 
所 降低 和 求解 时 间 的 有 效 减 少 ， 预 条 件 子 的 构造 往 
往 依赖 于 所 研究 的 物理 问题 [9], 对 预 条 件 子 的 改进 
是 提高 整个 数值 算法 性 能 的 关键 。 区 域 分 解 方 法 既 
可 以 作为 达 代 法 也 更 适合 于 作为 其 他 高 效 迭 代 求 解 
器 的 预 条 件 子 , 尤其 适合 于 当前 流行 的 分 布 式 并 行 
计算 机 环境 。 
基于 上 述 方法 ， 本 文 提 出 了 适用 于 求解 城市 大 
气 背 景 流 场 问题 的 Newton-Krylov-Schwarz (NKS) 
算法 。 使 用 非 精 确 Newton 算法 作为 非 线性 求解 器 ， 
在 每 个 非 精确 Newton 步 里 使 用 由 基于 区 域 分 解 方 
法 的 Schwarz 预 处 理 算 子 来 加 速 Krylov 算法 求解 
Jacobian 方程 。 其 详细 步骤 列举 如 下 : 
NKS 算法 ; 
步骤 1; 使 用 前 一 个 时 间 步 的 解 作 为 初始 值 
X' = X" 1， 
步 又 2: 对 = 01… 执 行 以 下 步 又， 直至 收敛 ， 
步 又 2. 1: 构造 Jacobian 矩阵 内 ， 
步骤 2. 2: 采用 带 有 预 处 理 的 Krylov 子 空 
间 方 法 (如 GMRES) 0 
了 (MD-IMXSY = —F"(XY), 
步 又 2.3: 通过 线性 搜索 获得 步 长 化 ， 
步骤 2. 4， 更 新 迭代 解 
XR = XR + TR 
在 此 步骤 中 ， 懈 为 待 求 非 线性 方程 组 芒 (X) 在 点 和 处 
的 全 Jacobian 矩阵 ，M 为 一 水 平 加 性 Schwarz 预 
条 件 算 子 ， 步 又 2. 2 中 非 精 确 是 指 求解 Jacobian 
II 
| TME) MS 十 了 (人 有 和 mm | F(X?) | 
式 中 水 为 线性 求解 器 的 相对 收敛 误差 。 
在 本 研究 中 ， 我 们 使 用 该 算法 求解 3.2 节 中 控 
制 方程 离散 后 形成 的 非 线 性 系统 (3)。 
本 章 介 绍 的 算法 基于 阿 贡 国家 实验 室 开发 的 开 
源 可 移植 并 行 软件 包 PETSc[10] 实 现 .PETSc 基于 基 
础 线性 代数 子 程序 库 (BLAS)， 线 性 代数 包 LAPACK 
以 及 消息 传递 接口 库 MPI 等 构成 线性 求解 器 、 非 线 
性 求解 器 以 及 时 间 步 进 积 分 器 三 个 主要 的 组 件 ， 为 
用 户 提 供 了 包括 求解 大 规模 稀疏 线性 方程 组 的 多 种 
Krylov 子 空间 方法 在 内 的 丰富 的 算法 以 及 预 条 件 
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计算 平台 上 具有 强大 的 偏 微分 方程 数 
值 求解 能 力 , 目 前 PETSc 已 经 成 功 应 用 于 优化 问题 、 
类、 计算 流体 动力 学 等 多 类 工程 与 科学 计算 
问题 。 本 文 的 数值 实验 国产 某 超 级 计算 机 上 展开 。 


5 数值 算 例 


作为 应 用 算 例 ， 我 们 对 真实 尺寸 的 深圳 地 王 大 厦 
附近 约 一 平方 公里 区 域 进行 高 雷诺 数 大 气流 场 进行 
计算 ， 通 过 流 场 信息 与 并 行 性 能 两 方面 对 算法 进行 
评价 。 本 节 采 用 NKS 算法 所 获取 的 结果 与 目前 主流 
的 商用 计算 流体 力学 软件 CFX 的 计算 结果 进行 对 比 
来 验证 算法 的 正确 性 。 为 使 模拟 结果 有 具有 可 比 性 且 
有 对 比价 值 ，NKS 算法 与 商业 软件 进行 数值 计算 时 
均 采 用 同样 的 几何 模型 、 流体 材 料 属性 、 计 算 区 域 、 
有 限 元 网 格 以 及 初始 与 边界 条 件 。 
5. 1 问题 描述 

深圳 市 红 岭 南路 、 宝 安南 路 路 、 滨 河 大 道 与 红 
宝 路 所 围 成 的 约 1000mx1000m 区 域 包括 多 座 深 
圳 地 标 建筑 , 如 地 王 大 厦 、 京 基 大 厦 以 及 万 象 城 等 。 
本 算 例 重 构 的 几何 模型 包含 了 该 区 域 大 部 分 的 高 层 
建筑 共计 29 座 ， 如 图 二 所 示 。 


> 


图 二 深圳 地 王 大 厦 附 近 区 域 建筑 模型 示意 图 


此 外 ， 在 剖 分 网 格 时 ， 为 保证 计算 精度 ， 生 成 
的 网 格 需 涵盖 所 有 建筑 群 的 主要 表面 特征 ， 即 所 谓 
的 包 面 技术 。 包 面 技术 在 进行 外 流 场 计 算 时 十 分 重 
要 ， 是 能 够 实现 复杂 模型 外 流 场 计 算 的 关键 所 在 。 
因此 ， 网 格 在 建筑 物 表 面 往往 需要 加 密 至 极 小 。 本 
文 将 包含 建筑 群 在 内 的 整体 计算 区 域 设 置 为 
Q= 4000mx4000mx800m ， 坐 标 原 点 位 于 京 基 大 厦 
与 地 王 大 厦 地 面 中 心 连 线 中 点 处 。 

非 结 构 四 面体 网 格 ™= 1 的 痢 分 是 进行 有 限 
元 计算 的 前 提 , 在 本 算 例 中 , 通过 商用 软件 ICEMCFD 
对 计算 区 域 Q 进行 非 结构 网 格 剖 分 ， 在 划分 网 格 
时 ， 我 们 对 建筑 物 壁 面 等 对 流体 运动 影响 较 大 区 域 
进行 加 密 处 理 。 经 剖 分 后 ， 本 算 例 的 非 四 面体 网 格 
总 数 约 为 1.3x10', 自由 度 约 为 DOF = 1.03 x 107， 
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建筑 物 壁 面 网 格 尺寸 达到 米 量 级 。 算 法 的 实现 需要 
对 计算 区 域 即 离散 后 的 网 格 进行 分 区 以 使 用 区 域 分 
解 方法 构造 适当 的 预 条 件 子 。 通 过 调用 开源 软件 包 
ParMETIS 进行 网 格 分 区 , 一 般 分 区 的 个 数 与 处 理 器 
个 数 相等 ， 且 各 个 分 区 所 计算 的 自由 度 个 数 近似 相 
同 以 达到 负载 平衡 ， 如 图 三 所 示 ， 不 同 颜色 的 网 格 
表示 在 不 同 处 理 器 上 计算 ， 需 要 注意 的 是 ， 本 示意 
图 并 非 将 真实 流体 计算 区 域 划分 而 是 将 建筑 物 作为 
固体 区 域 划 分 以 为 更 好 展示 并 行 分 区 概念 。 


图 三 计算 网 格 训 分 与 区 域 分 解 示 意图 

大 气流 体 取 25"c 标 准 大 气压 下 的 空气 ， 其 密度 

4198317 有 动力 粘度 系数 e934etyr。 设 i 
口 边界 风速 为 : 


Vs y V 
VW, = 一 (下 5) 太 十 (2 一 ae 
"(T8007 1800 


)1， 


其 中 设 V .=10m/s ， 时 间 步 长 AA=40s* ， 计 算 


1800s 内 流 场 的 变化 情况 。 算 法 中 的 线性 和 非 线性 
求解 器 的 相对 收敛 误差 分 别 取 为 1b* 和 10*。 分 别 取 


,10m1s 为 特征 速度 、 女 -50m 为 特征 长 度 ， 计 算 和 雷 


诺 数 为 : 
-| 
Re= purH _1.18Skg/m I -3236x107。 
HL 1.831x10 kg /ms 


作为 对 比 参 照 结果 ， 我 们 采用 商业 CFD 软件 
ANSYS. CFX 使 用 相同 网 格 数 、 边 界 条 件 的 问题 进行 
数值 模拟 ， 获 得 的 数值 结果 用 于 与 本 文 算法 获取 的 
数值 结果 进行 对 比 。 在 CFX 中 ， 瞬 态 项 采用 二 阶 向 
后 Euler 格式 , 均 方 根 残 差 (RMS) 收敛 误差 为 10 。 
5. 2 流 场 数值 结果 

图 四 、 图 五 所 示 为 1=1800s 时 刻 流 场 数值 结果 
对 比 。 从 建筑 物 表 面 压力 场 与 周围 速度 场 分 布 来 看 ， 
NKS 算法 获取 的 流 场 分 布 与 ANSYS 计算 结果 大 致 相 
同 ,在 建筑 物 迎 风 面 与 背风 面 分 别 测定 高 压 与 低压 ; 
在 建筑 物 之 间 的 街道 ， 空 气流 通 多 呈现 高 速 及 低速 
两 种 状态 ， 如 京 基 大 厦 附 近 空 气流 通 速度 较 高 ， 而 
不 同 ， 印 证 了 由 于 建筑 物 尤 其 是 高 层 建 筑 物 的 阻挡 
作用 ， 在 城市 建筑 群 密集 处 会 产生 “急流 ”或 “ 死 


流 ” 区 域 。 


图 四 压力 云图 对 比 〈 左 : ANSYS, 右 : NKS) 


图 五 Z = 80m 处 速度 场 分 布 云图 ( 左 : ANSYS, 右 : NKS) 


急流 区 空气 流动 剧烈 呈 汕 流 状态 ， 其 运动 状态 
十 分 复杂 ， 如 面 流 线 图 六 所 示 ， 为 空气 污染 扩散 模 
拟 等 进一步 分 析 带 来 了 一 定 的 困难 。 死 流 区 域 空气 
流速 低 ， 不 利于 空气 污染 物 的 扩散 ， 导 致 污染 物 聚 
集 形成 高 污染 区 域 。 因 此 ， 高 精度 的 大 气 环境 数值 
模拟 可 以 更 精确 的 获取 城市 某 一 小 区 域 的 空气 运动 
状况 ， 为 未 来 小 区 级 大 气 环境 与 质量 差别 化 预报 提 
供 了 技术 手段 。 


图 六 截面 Zz = 80m 处 面 流 线 图 


5. 3 并 行 可 扩展 性 结果 

本 小 节 给 出 基于 NKS 算法 求解 深圳 城区 大 气 环 
境 的 并 行 可 扩展 性 数值 结果 。Newton 非 线 性 迭代 与 
Krylov 子 空间 方法 线性 迭代 收敛 误差 分 别 设置 为 
10“ 和 10“， 数 值 结果 取 前 20 个 时 间 步 的 平均 值 。 


表 1 NKS 算法 的 并 行 性 能 


ILU LU 
p 

Newton GMRES Time Newton GMRES Time 
512 3.8 131.6 228.3 - - - 
1024 3.8 134.4 139.5 3.8 89.9 1051.7 
2048 3.8 143.3 93.2 3.8 114.2 351.2 
4096 3.8 149.2 66.1 3.8 141.3 117.3 
8192 3.8 154.2 42.7 3.8 168.2 60.1 


16384 3.8 168.6 31.8 3.8 204.8 37.6 


表 1 列 出 了 采用 ILU 和 LU 方法 作为 线性 求解 器 
时 ,不 同 处 理 器 条 件 下 计算 的 可 扩展 并 行 数值 结果 。 
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Newton，GMRES 和 Time 分 别 表示 为 每 个 时 间 步 内 非 [2] 《北京 城市 规划 建设 与 气象 条 件 及 大 和 气 污 染 


线性 迭代 次 数 、 线 性 迭代 次 数 以 及 计算 时 间 (单位 关系 研究 ) 课题 组 城市 规划 与 大 气 环境 . 北 
为 秒 )。 为 了 更 加 清晰 地 展示 可 扩展 性 性 能 ,本文 将 0 | 
放行 数值 结果 以 加 速 比 的 形式 列 于 图 8 中 。 泵 : 气象 出 版 社 , 2004. 


1051.7 32[ 


-1LU ET 
一 e 一 LU 16| 一 一 LU 


[3] T. Baileyd. User's guide for the industrial 
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data of dispersion from chemical warehouse 
fire. Atmos. Environ., 1999, 33(12): 
1937-1953. 


图 8 每 时 间 步 平均 计算 时 间 与 加 速 比 


从 并 行 测试 结果 我 们 可 以 看 到 ， 由 于 采用 了 前 
一 时 间 步 的 结果 作为 初始 值 ， 每 个 时 间 步 的 非 线性 
迭代 次 数 Newton 非常 少 , 非 线 性 迭代 次 数 几乎 与 处 [3] F. Chen. Developing an integrated Urban 


理 嚣 个 数 无 关 。 随 着 计算 核 数 的 增加 ， 线 性 迭代 次 modeling system for WRF. NCEP EMC 
数 GMRES 略微 增加 。 计 和 售 时 间 等 并 行 数值 结果 显示 Seminar, Camp Springs, 2005, 18. 


我 们 的 算法 具有 非常 好 的 可 扩展 性 ， 处 理 器 个 数 达 a 
6] 于 王 晓 云 直 模 卸 

0 到 万 核 时 并 行 效率 在 30% 左 右 , 而 采用 1 算法 作为 [G] 活 交 和 “大气 环 境 数 信 模 拟 在 城 

加 子 求解 器 时 可 获得 超 线性 的 加 速 比 。 市 小 区 规划 中 的 应 用 . 清华 大 学 学 报 (自然 科 


学 版 ), 2006, 46(9):1489-1494. 


结论 

Ny 0 所 化 [7] A..Martilll. Current research and future 
一 对 于 城市 大 气 环 境 模拟 问题 ， 高 保 真 度 往往 意 challenges in urban mesoscale modelling. Int. J. 
CA 。 味 着 庞大 的 计算 量 ， 几 何 多 尺度 、 高 雷诺 数 等 因素 Climatol., 2007, 27(14):1909-1918. 


致使 其 求解 非常 具有 挑战 性 ， 对 计算 机 的 硬件 与 算 [8] E. Batchvarova，9. E. Gryning. Advances in 


>< ie a Ane le 0 urban meteorology modelling. Advances in Air 
pe 法 以 获取 时 效 性 守 的 数值 结果 是 交 
= 值 模拟 0 关键 i 基于 国产 大 规模 计算 未 统 > Pollution Modeling for Environmental Security. 
过 一 Lo 人 ~ 二 口 ， 
r= 本 文通 过 对 大 规模 非 线性 系统 进行 高 效 求解 器 和 预 Springer Netherlands, 2005: 23-32. 
处 理 技术 的 研究 ， 提 出 一 种 可 扩展 并 行 [9] Knoll D and Keyes DE. Jacobian-free 
Newton-Krylov-Schwarz 算 法 。 非 线性 方程 采用 非 精 Newton-Krylov methods: A survey of 
确 Newton 方 法 进行 求解 , 在 每 个 Newton 步 , Jacobian 加 
系统 通过 基于 区 域 分 解 方法 的 限制 加 性 Schwarz 预 approaches and applications [J]. J. Comput. 
条 件 子 处 理 ,然后 使 用 以 GMRES 方 法 为 代表 的 Krylov ee 
子 空 间 迭 代 法 作为 线性 求解 器 进行 求解 。 作 为 应 用 ， [10] Balay S，Buschelman K，Gropp WD, et al. 
我 们 对 深圳 地 王 大 厦 附 近 约 一 平方 公里 区 域 进行 大 PETSc Users Manual [R]. Argonne National 
气流 场 进 行 计算 。 数 值 结果 显示 ， 本 文 的 算法 在 扩 Laboratory, 2016. 


展 至 万 核 处 理 器 平台 时 仍 具 有 非常 良好 的 可 扩展 并 
行 性 能 ， 为 未 来 对 整个 城市 区 域 进 行 高 精度 如 网 
格 分 辨 率 达 到 米 级 ， 未 知 数 个 数 达 十 亿 量 级 ) 的 大 
气 环境 模拟 提供 了 一 种 可 参考 的 算法 。 
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